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Abstract 



A generalized stochastic method for projecting out the ground state of the 
quantum many-body Schrodinger equation on curved manifolds is introduced. 
This random-walk method is of wide applicability to any second order differ- 
ential equation (first order in time), in any spatial dimension. The technique 
reduces to determining the proper "quantum corrections" for the Euclidean 
short-time propagator that is used to build up their path-integral Monte Carlo 
solutions. For particles with Fermi statistics the "Fixed-Phase" constraint 
(which amounts to fixing the phase of the many-body state) allows one to 
obtain stable, albeit approximate, solutions with a variational property. We 
illustrate the method by applying it to the problem of an electron moving on 
the surface of a sphere in the presence of a Dirac magnetic monopole. 
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I. INTRODUCTION 



The importance and difficulty of solving models of interacting quantum particles is hard 
to overstate. It is well known that the correlated motion of those particles gives rise to a 
wide variety of physical phenomena at different length and time scales, spanning disciplines 
like chemistry, condensed matter, nuclear, and high energy physics. Novel complex struc- 
tures can emerge as a consequence of the competing multiple-length scales in the problem. 
Nonetheless, only a reduced set of interacting problems admits exact closed form solutions 
HI] and the use of numerical techniques becomes essential if one is looking for accurate solu- 
tions not subjected to uncontrolled approximations. Among those techniques, the statistical 
methods 0] offer the potential to study systems with large number of degrees of freedom, 
reducing the computational complexity from exponential to polynomial growth. This scaling 
behavior is particularly relevant when one recognizes that most of the interesting phenomena 
in many-body physics occurs in the thermodynamic limit Unfortunately, for fermions 
(i.e. quantum particles obeying Fermi statistics) the sign problem plagues all useful stochas- 
tic algorithms and causes the variance of computed results to increase exponentially with 
increasing number of fermions [Q. 

On the other hand, the growing interest in physical systems whose state functions are 
defined on a general metric space makes the quantum mechanics of interacting particles 
in curved manifolds no longer a mere intellectual exercise, but one with very practical 
consequences. Perhaps the most well-known examples can be found in cosmology (e.g., 
matter in strong gravitational fields, atomic spectroscopy as probe of space-time curvature 
1^), but the subject is certainly not exclusive to this field. In condensed matter a very 
elementary case is provided by a deformed crystal. Less well-known ones are mesoscopic 
graphitic microtubules and fuUerenes. All these physical systems are ubiquitous in nature 
and the crucial role the curvature of the manifold plays has been confirmed by experimental 
observations (e.g. spectrum of collective excitations 0). Therefore, the development of 
stable quantum methods with polynomial complexity in Riemannian manifolds represents a 
real challenge for many-body theorists. 

The present manuscript deals with the (non-relativistic) many-particle Schrodinger equa- 
tion in a general metric space and its solution using stochastic techniques. In particular, we 
will show how to construct approximate solutions (wave functions) for systems with broken 
time-reversal symmetry (e.g. electrons in the presence of external electromagnetic sources) 
avoiding the infamous "phase problem" 0. The main difficulty is to define a probability 
measure (semi-positive definite) which allows one to get the complex-valued state with no 
asymptotic signal-to-noise ratio decay in Euclidean time. This translates into a problem of 
geometric complexity, which is solved approximately using constraints in the manifold where 
the wave function has its support. In this way, we get stable but approximate solutions which 
can be systematically improved. 

Among the large variety of problems one can attack, we decided to choose the general 
problem of fermions in the presence of external gauge fields to illustrate the main ideas. 
The effects of an external magnetic field on a system of electrons can be profound 0]. The 
field couples to the electron's charge and spin, modifying its spatial motion and lifting its 
spin degeneracies. The field can also create spatial anisotropy, effectively reducing the di- 
mensionality of the system from three to two. The combination of the reduced dimension 
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and the field itself is known to have novel consequences. For example, in a system of non- 
interacting electrons hopping on a square lattice, the field transforms the energy spectrum 
from the simplicity of trigonometric functions to the complexity of a field-dependent self- 
similar structure (Hofstadter's butterfly) whose depth mathematicians are still fathoming 
0]. The combination of the reduced dimensionality, strong particle interactions and the field 
itself is known to have novel consequences, like the formation of isotropic fractional quan- 
tum Hall fluids |T^, which are incompressible states of the two-dimensional homogeneous 



Coulomb gas. 

The projector (zero temperature) method we will introduce uses random-walks to solve 
a general multidimensional partial differential equation second order in space coordinates 
and first order in time. Whenever mention is made of a random-walk we mean a Markov 
chain that is defined as a sequence T^i, 7^.2, ■ ■ ■ , TZr of K random variables that take values 
in configuration space, i.e. the space of particle positions. As usual, what characterizes a 
random- walk is its initial probability distribution and a conditional probability that dictates 
the transition from TZi to TZi+i- This transition probability is non- unique and discretization 
dependent Among all the possible choices we will require a prepoint discretization of 
the transition probability (short-time propagator) because we will use Monte Carlo methods 
to generate the walkers. 

The paper is organized as follows. In Section |T| we present the formulation of the 
general problem of fermions in curved manifolds. In particular, for illustration purposes 
and to fix notation, we develop the formalism for spin-| particles in the presence of an 
external electromagnetic potential. Then, we show how to project out the lowest energy 
state of a given symmetry in a manifold with curvature, and discuss the resulting Fokker- 
Planck equations for various distribution functions. Once the problem is precisely defined 
we develop, in Section |T|, path- integral solutions to those multidimensional differential 
equations, and give an interpretation of the emergent "quantum corrections" in the Euclidean 
action. The path-integral solutions are evaluated using Monte Carlo techniques in Section 



IV. There, we provide an stable step by step practical algorithm which emphasizes the 



subtle changes (with respect to the standard Diffusion Monte Carlo (DMC) technique) due 
to the metric of the manifold. In Section we apply such computational implementation 
to the problem of an electron moving on the surface of a sphere in the presence of a Dirac 
monopole. Finally, Section |VI| summarizes the main findings and discusses the relevance of 
the stochastic method as applied to the physics of quantum Hall fluids. 



II. FERMIONS ON RIEMANNIAN MANIFOLDS 

Notation. Consider a differentiable manifold Ai of dimension d (e.g., for the two-sphere 
S^, d=2) with coordinates = {x}, . . . ,xf) defined on it. If is a Riemannian manifold, 
then it is a metric space, with metric tensor g^'^{ri) = g'^'^{i), such that the distance ds 
between two points in Ai is ds'^ = gp.v{i) dx^dx'^ in the usual way ||12|. The metric tensor 



is positive definite and symmetric g^'^ = g'^^ (as we will see, this condition is important to 
define a probability density distribution), and is a function of the coordinates with the 
property g^-^g'^^ = 5^. Let us consider the coordinate transformation h: xf = h^{x'l, . . . , x'f). 
Then, a generic second order contravariant {T^'^) and covariant tensor {T^^) transform as 
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rp^V _ ^-^i ^-^i rplaP rp _ ^-^i ^-^j rpl / -,\ 

dxTdxf '^^ dxt dx\ '^^ ' ^ ' 

respectively. Throughout the paper Einstein's summation convention for repeated indices is 
assumed (/i, z/ = 1, . . . , ci). 

Formulation of the problem. In this article we will be concerned with finite interacting 
fermion systems in the presence of an external electromagnetic potential a^(rj) = a^{i) = 
(A(i), = 0) (B = V A A represents a uniform field, A and are the vector and scalar 
potentials, respectively) whose Hamiltonian for motion on the manifold, in the coordinate 
representation, is given by 

H = eo + r({r,},{s,}) (2) 

with 

ph ^ p2 AT 

Ho = -DA + E [2«"(05m + 9-'^'(^)d, {g'/'(^)a'^{^))] + E , (3) 



i=l ^"^ i=l 



2m*c 

(9^ = d/dxi, and A = X^^Ii ^{i)-, where 

A = g-"% (g^'^g'/'d^ (4) 



is the covariant Laplace-Beltrami operator and V is a potential energy operator. Notice 
that we use the conventional notation where the transformation between different forms 
of a given tensor is achieved by using the metric tensor (e.g., a'^ = g^^a^, = gf^i^a''), 
and g^^'^ = Jdet g^u- This Hamiltonian characterizes the dynamics of non-relativistic 



indistinguishable particles of mass m*, charge e and spin Sj = | in a curved space with 
metric tensor g^^ , and D = h^ /2m*. We have assumed that the quantum Hamiltonian IH 
in curved space has the same form as in flat space (this amounts to a particular operator 
ordering prescript ion . ) 

Given the previous ordering, one can rewrite the Hamiltonian above in terms of the 
generalized (hermitian) canonical momentum = —ih(d^ + |9^(ln (7^/^)) 

_ 1 ^ 

H = ^ E 9~'^\^) n, g'/\z) g^'^m^ g-''\{) + V ({r,}, {.,}) , (5) 

where the kinetic momentum 11^ = P/^ — f O/x- The first term in Eq. ^ represents the kinetic 
energy of the system and is the non-relativistic approximation to the Dirac operator. V 
includes the sum of one and two-body local interaction terms (and background potential in 
the case of a charge neutral system) and Zeeman contribution. The potentials are assumed 
to be finite almost everywhere and can only be singular at coincident points (r^ = rj, V i 7^ j) 
We are interested in the stationary solutions of the resulting multidimensional Schro- 
dinger equation 

ih dt\^) = e 1^) , (6) 

and will restrict ourselves to Hamiltonians IH which are time-translation invariant. In the 
usual space-spin formalism the A^- fermion states characterizing the system, {X\^) = ^^(X), 
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and all its first derivatives belong to the Hilbert space of antisymmetric (with respect to iden- 
tical particle (rj, Sj)-exchanges) square- integrable functions Hn = C^iM^) ® C^^, defined 



as 



n 



N 



1^ I p-.xif = , and ll^ll = vW*) < 00} , (7) 



where X = (JZ, S) {Ti = (ri, . . . , fat) and S = (cxi, . . . , gn) are discrete spin variables) and 
Pij represents the permutation of the pairs (r^, cTj) and (rj, aj). is the Cartesian product 
manifold of dimension dN. 

Since the system Hamiltonian can be written as IH = ]H.ti{TI) + ]H.s{^)j the last term 
representing the Zeeman coupling, the many-body wave function \E'(7^, S) can be expressed 
as a tensor product of a coordinate and a spin function (or a linear combination of such 
products), 

^(7^, r) = $(7^) ® e{s) . (8) 

We want to construct A^-fermion eigenstates of H, \I' , that are also eigenfunctions of the 
total spin 5"^ (5* = Y^iLi Si), 

^(X) = s{s + 1) ^(X) , (9) 

and this is always possible since H, 5*^ = 0. Thus, the configuration part ^(Jt) must 
have the right symmetry in order to account for the Pauli principle. It turns out that a 
coordinate state $(ri, . . . , r^, r^+i, . . . , r^r) which is symmetrized according to the Young 



scheme ||13[ and has total spin s = y — will be antisymmetric in the variables ri, . . . , r^. 




and antisymmetric in the variables r^+i, . . . , tat. Moreover, $ possesses the property of 
Fock's cyclic symmetry, 

/ N \ 

= , (10) 

where, in this case, Pkj refers to the transposition of particle coordinates and r^. This 
last condition is a very useful one for testing the symmetry of a given coordinate function. 

Quantum projection on curved manifolds. For a given total spin s we are thus left with 
the task of solving the stationary many-body Schrodinger equation ]H.n^{TZ) = E^[7V), 
where $(7^) = (7?.|$) satisfies the symmetry constraint discussed above. In particular, we 
are interested in the zero temperature properties of this quantum system, i.e. its ground 
state properties. To this end, we study the Euclidean time evolution of the state i.e. we 
analytically continue Eq. |^ to imaginary time (Wick rotation, t — —ith) 



$ , (11) 



whose formal solution $(t) = h({t) $t = exp[— t(lll7^ — Et)]^t is used to determine the 
limiting distribution 

$0 oc /im <l>(t) , (12) 
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which is the largest eigenvalue solution of the evolution operator U{t) compatible with the 
condition (<I>o|'^'t) 7^ 0, where $t is a parent state and Et is a suitable (constant) energy 
that shifts the zero of the spectrum of IH7J. 

We would like to solve the multidimensional differential equation Eq. |TT| using initial 
value random walks. In this way, starting with an initial population of walkers (whose state 
space is A4^) distributed according to piJZ, t = 0) = $1- ($7^ must be positive semi-definite), 
the ensemble is evolved by successive applications of the short (imaginary) time propagator 
U{t) [t = t/M, and M is the number of time slices) to obtain the limiting distribution $o- 
Then, we can introduce a "pseudo partition function" 

Z = U{t) $r) (13) 

in terms of which we can determine the ground state energy Eq as 

Eq-Et= lim -i InZ . (14) 

Similarly, other ground state expectation values, e.g., {^o\0 $0)) can be obtained as deriva- 
tives (with respect to a coupling constant J) of a modified pseudo partition function Zj 
whose evolution operator has a modified Hamiltonian, IHt^ + JO. 

In order to reduce statistical fluctuations in the measured quantities (i.e., observables) 
one can guide the random walk with an approximate wave function, $g, which contains 
as much of the essential physics as possible (including cusp conditions at possible singular- 
ities of the potential V). Then, instead of sampling the wave function <l>(t) one samples 
the distribution f{TZ,t) = ^{t)^G (properly normalized) with the initial time condition 
f{7l,t = 0) = ^T^G- Expectation values of operators O (observables) that commute with 
the Hamiltonian have a particularly simple form for guided walkers. For instance, 

hm ^ ($5^a$.),,_, , (15) 

where the average {A) j stands for 

/^.a;/(7^,t^oo) ' 

/(??., t 00) is the long-time stationary probability of the system, and the (invariant) 
volume element u: is given by the ciA^-form 



■ N 
.1=1 



dxi A ■ ■ ■ dxf A ■ ■ ■ dx^ . (17) 



Remember that in a general metric space the resolution of the identity operator with respect 
to the spectral family of the position operator is 

i=[ u:\n){n\. (18) 

JM" 
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It is important to stress that $r and the guiding function $g can, in principle, be different 
functions, although most of the practical calculations use the same function. It turns out 
that this importance sampling procedure is decisive to get sensible results when the potential 
V presents some singularities. 

Notice, however, that the quantum Hamiltonian H breaks explicitly time-reversal sym- 
metry, meaning that in general $ will be a complex-valued function. Even if $ were real- 
valued, because it represents a fermion wave function it can acquire positive and negative 
values (the case where H(Z') is totally antisymmetric being the exception). Then, it is clear 
that we cannot in principle interpret $ or / as a probability density. 

For reasons that will become clear later we will be interested in sampling the probability 
density f{lZ,t) = \f(7l,t)\. The generalized diffusion equation in curved space for the 
importance-sampled function / can be derived directly from Eq. ^ with the result 

_ N 

dj =DY: [9-'^^^)^, {g^''{i)g"'{i){dj-f F,))] - {El - Er)f , (19) 

i=l 

where the drift velocity Fi,{TZ) = duln^Q, and the "local energy" of the effective ("Fixed- 
Phase") Hamiltonian Hpp (see Eq. ^ and its derivation) is El{TZ) = Hfp^g with 



^ r e e 

Hpp = -DA + DY: (d^xin) - -a>^(t))(d,xin) - -a,{i)) 



+ v{n) , (20) 



where x('^) is the phase of the many-body state i.e. $ = |$|exp[ix]. The differ- 
ential equation satisfied by the distribution function / is formally equivalent to the one 
describing Brownian motion on a general manifold (including generation and recombination 
processes), and corresponds to a Kramers-Moyal expansion with exactly two terms. In fact, 
we can rewrite the equation above as a Fokker-Planck equation for dN continuous stochastic 
variables {rj}j=i^...^7v 

dtf = [C^^-{EL-ET)}f, (21) 
where the (time-independent) Fokker-Planck operator C-p-p is given by 

TV 

£fp« = E [9,du {D'-'ii}*) - d, (D'^i'R)*)] . (22) 

1=1 

The diffusion matrix (contravariant tensor) D^'^ and drift Z)^ (which does not transform as 
a contravariant vector) are given by 

D^"" = D g^"" (23) 
= D^^F^ + d^D^" - D^'^'Tl^ , (24) 

where P^^^ is the Christoffel symbol of the second kind 

r^. = \g'"' {d.gup + d^g^p - d,g,,) (25) 
T:^ = ^d,\ng, (26) 
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and the modified local energy 



(27) 



Notice, however, that singularities in the "quantum corrections" [|14[ to the local energy 
El due to the metric, can induce very large fluctuations in E^. Moreover, the probability 
density / does not transform as a scalar function (/(7^, t) uj = f{TZ', t) u}', where the primes 
represent the transformed coordinates and lj = dx\ A ■ ■ ■ dxf A ■ ■ ■ dx% is a volume element 
in Ai'^). Therefore, it is more convenient to work with a probability density that is a scalar 
and which is defined as 



f{n,t) 



N 



U9 

.i=l 



1/2/ 



m,t) 



(28) 



The differential equation / satisfies is of the form Eq. ^ with bar quantities replaced by 
unbar ones (e.g. £fp — > Cpp). It turns out that D^'^ = D^'^ and the drift (which is not a 
tensor) 



(29) 



Note that in this case the quantum correction to the local energy vanishes. Furthermore, 
if the metric is diagonal, i.e. g^^ = g^^'^5^u-, then the correction to the fiat space drift 
also vanishes, i.e. d^D^^^ + D^^^V^^ = 0, and D'^ = Dg~^/'^F^. This last remark is quite 
important, specially for = 2 where it is always possible to choose a coordinate system 
(rj = where the metric tensor is diagonal (conformal gauge [0), and use the 

conformal parameterization {zi = Q + i^f, zi = ^j—i^f ) which greatly simplifies the resulting 
expressions (see Section 0). 



III. PATH INTEGRAL SOLUTIONS 



The generalized Fokker-Planck Eq. ^ describes the time evolution of a distribution 
function / which is completely determined by the distribution function at t = to = 0. 
In this sense it describes a continuous stochastic process that is Markovian. Because it 
represents a Markov process, the conditional probability that if the system is in TZ at time 
t = it will jump to TZ' in time t (importance-sampled Green's function) G{TZ — > TZ'; t) 
contains the complete information about the process, and it follows that the probability 
densities fiJZ^t + r) and f(7l,t) are connected by 



/(7^', t + T)= f ujG{n^ 7^'; r) /(7^, t) 



(30) 



where the Green's function G{TZ TZ'',t) is a transition probability for moving particles 
from 7^ to 7^' in time r with the initial value G(7^ W; 0) = [uf=i 9~^^^{i)] S{n - TZ'), 
and is formally given by 



G{n ^ 7^'; r) = <I>G(7^') (7^'| exp[-T{HFp - Et^ C^) 



(31) 



with 
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N 



$G(7^) 



.1=1 



1/2/ 



$G(7^) 



(32) 



Path integral solutions for / may be derived from the transition probability density for 
small T. Iteration of the Chapman- Kolmogorov equation for G allows one to express the 
evolution of f{TZ',t) from the initial distribution f{TZ,t = 0) in terms of the short-time 
Green's function as 



[ cvo G{TZm-i ^ 7^M; r) ■ • ■ G(7^o ^ 7^l; r) /(7^o, 0) , (33) 



where t = Mr; we identify TZq = TZ and TZm = TV ■ By simple inspection the solution of 
the generalized Fokker-Planck equation / stays positive if it was initially positive (i.e., if 

/(7^,o)>o.) 

It is clear then that the functional integral representation of / requires knowledge of the 
infinitesimal evolution operator. It is well-known in a similar context |Tl]] that the integrand 
of the functional integral is not unique, it is discretization dependent (compatible with the 
Markovian property of the paths since 'R-{t) is only sampled at t — r and t.) On the other 
hand, it is crucial for numerically simulating those paths to use a discretization where the 
drift velocity and diffusion are evaluated at the prepoint in the integral equation. Following 
Feynman |jT6| we have determined the functional form of G{TZ —>■ TZ'; r) to 0(t^). We are 
going to present the final result and omit the details of the calculation which are just sim- 
ple (although lengthy) manipulations of Taylor series expansions and gaussian integration. 
Thus, to O(r^) the short-time conditional probability is given by 



N 



G{TZ 7^'; r) = G'^,(7^ ^ 7^'; r) \{ G^.iJZ 7^'; r) 

i=l 



where 



G^,(7^ ^ 7^'; r) = exp 



(34) 



(35) 



and 

G°(7^^7^';r) 



AnDi 



d/2 



exp 



X. 



ADt 



, (36) 



that is, a gaussian distribution with variance matrix 2D^^ and mean + tD^(JZ). Sample 
trajectories (continuous but nowhere different iable) are generated by using the Langevin 
equation associated with the process, i.e. x^" = xf + tD^{1Z) + y/r rj , where 77 is a gaussian 
random variable with zero mean. Note that D'^ and D^^ are evaluated at the prepoint in 
the integral equation. 

Therefore, the Wick-rotated path integral for G{1Z — * TZ'] t) is 



Gijl 7^'; t) 



n{t)=n' 
n{o)=n 



(37) 



where the measure V[(jj{t)] = limM^oo(4vrDr) ^^"^/^ uji ■ ■ -ujm-i, and the Euclidean action 
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N 

S [7^(t)] = i dt' \ ^ ^(xf - D^m')]) g,^{i) {x\ - D^\nm + E^\n{i!)\ - Er 




(38) 



The integrand above represents a generalized Onsager-Machlup function |T^, and the dot 
is a short-hand for time derivatives. In closing this Section, we would hke to mention that 
functional integral solutions for / can be obtained from previous expressions after making 
the replacement {D^.E^) [D^.Ej^). 

Interpretation of the quantum corrections. The general methodology we have developed 
so far, can be equally applied to other situations which do not necessarily involve a curved 
manifold such as, for instance, particles moving in a medium with position dependent diffu- 
sion constant. In order to adapt our previous formalism, we need to understand qualitatively 
the origin of the quantum corrections to the short-time propagator obtained above. To this 
end, we will illustrate the general idea with the following Id equation {A4 = IR, = 1) 

dliDix)f{x,t)) = dtfix,t) . (39) 

The "standard" approach to finding the Green's function for this problem is simply to solve 
the equation dl{D{x)G{x,t)) = dtG{x,t) subject to the boundary condition G{x, 0) = S{x). 
We can do this by taking the Fourier transform in x: 

- — / dk' D{k - k') G{k\ t) = dtG{k, t) , (40) 

with the boundary condition G{k,0) = 1. This is more complex than the usual diffusion 
equation because we get a convolution of D{k) and G{k,t). However, we can find an ap- 
proximate solution for G{k,t), valid for small time, by noting that the boundary condition 
implies that, for small times t, G{k,t) — G{k',t) ~ 0{t) and so we have 



/ dk' D{k - k') Gik, t) + 0(t) = dtG{k, t) (41) 

27r JM 



k' 

IM 

or, simply, —k^D{x = 0)G{k,t) + Oit) = dtG{k,t), where D{x = 0) is D{x) evaluated at the 
prepoint. For small times we can ignore the order t term and solve for G with the result 
G{k,t) = exp[— t k'^D{0)]. Taking the inverse Fourier transform we finally have 

G{x, t) = , ^ exp\-xy4D{0)t] + 0{t^) , (42) 
^J4nD{0)t 

which is just the plain Green's function for a Id random- walk with D{x) evaluated at the 
prepoint and with no quantum corrections. This is, in fact, the result that the Green's 
function for the Jacobian times / has no quantum corrections. 

To make things clear let us look at a different equation with the same boundary conditions 

D{x)dlG{x, t) = dtG{x, t) , (43) 

which characterizes a free Brownian particle in Id. Again, taking the Fourier transform we 
get 
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- -5- / dk' k'^ D{k - k') G{k', t) = dtG{k, t) . (44) 
zvr JM 

Making the same approximation as above, replacing G{k',t) with G{k,t) in the integrand, 
making an error of 0{t), and noticing that 

^ f dk' {k - k'f D{k') = eD{0) + 2ki d^D{0) - dlD{0) (45) 

we get G{k,t) = exp[— t {k'^D{0) + 2ki dxD{0) — dlD{0))]. When we take the inverse Fourier 
transform the terms with derivatives of D{x) ai x = give precisely the quantum corrections 
for this simple case 

G{x, t) = , ^ exp[-(x - 2tdxD{<d)f/AD{<d)t + t 9^D(0)] + 0{f) . (46) 
/47rL'(0)t 



This is the result that the Green's function for / has quantum corrections. 

The general case is just as simple, and the following rule emerges. Given any second 
order differential equation (first order in t), no matter how many dimensions, with or without 
curvature, with or without a position dependent diffusion constant, the rule for obtaining 
the short-time Green's function with everything evaluated at the prepoint is as follows: 
Bring all derivatives in each term all the way to the left of that term. Once this 
is done one can simply write down the Green's function for a generalized diffusion process 
assuming the D{x), D^, D'^'^ , whatever position dependent terms they may be, are constant 
and evaluated at the prepoint. The quantum corrections are then seen to be simply those 
extra terms we get when commuting the derivatives to the left. 

Fermion-phase problem: Fixed-Phase method. It is evident that one cannot make a 
probability density out of a complex and/or antisymmetric wave function. This is the 
reason why we decided to write down Fokker-Planck equations for the distribution / (or /) 
and not /. Nevertheless, the phase factor associated with the original complex distribution 
must show up in the evaluation of the expectation values. It is well-known that this causes 
the variance of the computed results to increase exponentially with increasing number of 
degrees of freedom. This problem is known as the fermion-phase |]TB[ catastrophe, and 



in this Section we will review a method ^ to obtain stable, albeit approximate, path- 
integral solutions whose stochastic determination has a polynomial, instead of exponential, 
complexity. The generalization of the ideas presented below to bosons (with complex- valued 
states) or anyons in general is straightforward |T9(| . 

In order to avoid cumbersome notation which would obscure the main ideas, here we will 
consider a simplified version of our original Hamiltonian 

^^ = ^ + ^(^) ' (47) 

where, for simplicity, we introduced the vector notation, 11 = (11(1), ■ ■ ■ , n(A^)), and in 
this subsection the same convention will be used for other bold quantities. The microscopic 
equations governing the imaginary-time evolution of our interacting system can be found 
from a variational principle of the form 6S[TZ(t)] = 0, where the Euclidean action is given 
by 
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^[7^(t)] = ^*dt'^^cj|^[n<i>]*- [n$] + t/ $*<i> + <i>*9t/<i>| . (48) 

Finding the states $ and $* which minimize S[TZ{t)] is equivalent to solving the 
Schrodinger equation and its complex conjugate. Equivalently, we can consider as inde- 
pendent real fields the phase and modulus of the wave function, i.e. |$| and x such that 
$ = |$|exp[ix], and perform independent functional variations on S[JZ(t)]. The resulting 
Euler-Lagrange equations are: 

-dt\<!>\ = ^e{exp[-tx]Mn(^} = Hfp 

(49) 

^ -(dtx) = Qm[exp[-tx]Mn<^}\<^\ = -| d^ J'' 

where 3?e and 53m stand for the real and imaginary parts of the expressions in brackets, 
respectively, 

= 1^ + ^ = ^ + ^ (^'^(^) - ^«') ■ (^MX(7^) - ^a,) , (50) 



and 



2m* 



^ ($* [n^<i.] + [n^<i.]* $) = ^ (nd.xin) - ^a^) (51) 



is the probability current. The singularities of x(^) occur at the zeros of $, which generically 
have CO dimension two. 

So far we have simply mapped the original fermion problem into a bosonic one for |<l>| 
but still coupled to its phase fiuctuations. Alternatively, one can regard this as a gauge 
transformation of the original fermion problem, whose effect is to add a non-local gauge 
field potential, d'^X: giving rise to a fictitious magnetic field. Notice that it is this gauge 
field that contains information on particle statistics. Moreover, although the geometry of 
x(7^) can be altered by a gauge transformation, the singularities remain invariant. 

The Fixed-Phase (FP) method [0] consists in making a choice for the phase, xt, and 
solving the bosonic problem for |<l>| exactly using stochastic methods. The method is stable 
and has the property of providing a variational upper bound to the exact ground energy 
Eq, Epp > Eq (the equality holds when xt is the exact ground state phase), and for a 
given Xt the lowest energy consistent with this phase. The trial phases xt should conserve 
the symmetries of H-;^ and particle statistics (for time-reversal invariant systems there is 
a way for systematically improving a given mean-field phase using projection techniques 
pop. Notice that the FP method projects out the lowest energy state of a given symmetry. 
Therefore, the method allows one to compute also excitations which are "ground states" of a 
particular symmetry. For ground state properties of real symmetric Hamiltonian operators 
the FP approach reduces M to the Fixed-Node method [ETI]. 



IV. COMPUTATIONAL IMPLEMENTATION 

In this Section we present an algorithm for computing the ground state properties of 
quantum many-body systems defined on a curved manifold with general metric g'^'^. As 
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mentioned in the Introduction, this can be accomphshed by performing all multidimensional 
integrals using Monte Carlo techniques. In this way, ground state expectation values are 
obtained by averaging over a large number of particle configurations generated according to 
a certain limiting probability distribution p{TZ, t oo). There is some freedom in the choice 
of this distribution p(7l,t), however, to reduce statistical fluctuations in the observables to 
be computed it is more efficient to use the so-called importance- sampled distribution / (JZ, t) , 
which is the product of the absolute value of the solution of the time dependent Schrodinger 
equation $(7^, t) and some positive function $g(7^) that is the best available approximation 
to the modulus of the ground state eigenfunction. In a curved manifold, on the other hand, 
it is more convenient to work with the modified importance- sampled distribution f{Tl,t), 
which is defined as a product of the conventional importance-sampled distribution / (JZ, t) 
and the metric (see Eq. |28| ). 

The propagation of particle configurations in time r is determined by the conditional 
probability (Green's function) G{1Z 1Z']t), whose separation into a diffusion (plus drift) 
and branching parts (see Eq. |3^) makes it very simple to simulate numerically. The gaussian 
term represents propagation according to the equation x'-^ = + rD'^ilZ) + ^/T rj , where 
?7 is a gaussian random variable with zero mean. The effect of the term tD'^{TZ) is to 
superimpose a drift velocity on the random diffusion process so that particle configurations 
are directed towards regions of configuration space where $g('^) is large. The branching 
term Gb{TZ — > TZ';t) in Eq. determines the creation and annihilation of configurations 
(walkers) at the point TZ' after a move. If the size of the ensemble of walkers at any time t 
is defined as 

then, its rate of change is given by 

dtV{t) = - [ uj [EL{n) - Et] f{n, t) . (53) 

Therefore, if the local energy El{TZ) is a smooth function of TZ, and the trial energy Ej^ is 
suitably adjusted, the size of the ensemble of walkers will remain approximately constant as 
the configurations propagate. In particular, if the local energy is constant and equal to Et 
then the fluctuations in the ensemble size will vanish. To ease notation, in the rest of the 
paper we will only consider the standard situation $t = $g- In such a case, ground state 
expectation values of a generic observable O will be computed as 

(<I>g| W(t) <I>g) c-//!*— ) J^^ V{t^oo) ^ ^ "-^^ ^ ^ ' 

In the following we present a step by step computational algorithm for implementing the 
stochastic approach discussed above. Most parts of the algorithm follow closely the standard 
DMC method, described for instance in |21[] , but we believe that it is useful to present these 
steps in detail here because a number of straightforward but subtle modifications due to the 
space curvature are involved. To be more specific (without loosing the general features), let 
us present the algorithm as it is applied to fermionic systems within the FP approach. 

Algorithm. 
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§ 1. Construct a guiding wave function $g, which is the best available approximation to 
the exact ground state (or lowest energy state of a given symmetry). In principle any choice 
of $G which has finite overlap with the exact state is acceptable, but the more accurate 
is the faster the convergence to the stationary solution will be, and the lower the statistical 
fluctuations will be as well. Recall that zero variance is obtained when the guiding wave 
function is equal to the desired ground state (and the ground state is bosonic). 

§ 2. Given the guiding function compute the quantum drift velocity Fi, = d^ln^Q and 
the local energy = ^q^Hfp^g^ where Hpp is the FP Hamiltonian defined in Eq. [2^. 
is used to construct the drift vector according to Eq. and the drift and local energy 
are used in evaluating the short-time Green's function according to Eqs. and |36 . 

When a simple guiding function can be constructed, the expressions for the drift and the 
local energy can be evaluated analytically, but in most cases this must be done numerically. 

§ 3. A set of initial configurations or walkers {TZjit = 0)} (j = 1, ■ ■ ■ , iV^) is created, 
such that particles in each walker are distributed according to the modified importance- 



sampled distribution f{TZ,t = 0) 



Variational Monte Carlo (VMC) technique). 

§ 4. Each walker TZj is diffused for a time r according to the gaussian part of the prop- 
agator Hill Gl{lZj — > IZ'j] t). This can be accomplished by moving each particle coordinate 

according to 



(which is equivalent to using the 



X. 



Xi 



(55) 



where ?7 is a Gaussian random variable with a mean of zero and a variance of 2D'^'^ = 2Dg^" . 
§ 5. The move from TZj to TZj is then accepted with a probability 



A{TZ, 7^' ; r) = min(l, W{TZ,, TZ'- r)) 



(56) 



where 



W{TZ,TZ';t) = 



N 



$G(7^' 



<fG(7^) 



GjTZ' TZ- t) 
G{TZ TZ'; r) 



(57) 



This step ensures detailed balance in the Monte Carlo procedure. A typical acceptance ratio 
is in excess of 99%. Notice that if G{TZ TZ'; r) is the exact Green's function, and not its 
short-time approximation, then W is unity and this step is not necessary. This is because the 
Green's function of an Hermitian operator is symmetric, i.e. G{TZ — > TZ'; r) = G{TZ' ^ TZ;t), 
but this is not the case for any importance-sampled distribution function equation. 

§ 6. After all the particles of a given walker have been diffused from the initial position 
TZj to the position TZ'j and the move is accepted, then the values of the drift D'^, the 
local energy El, and $g are updated. We could have equally well move all particles at 
once in step 4 before step 5, however, the acceptance probability for a given time step is 
reduced considerably. Depending upon the particular problem one can adopt one of the two 
strategies: single or multiparticle moves. 

§ 7. The multiplicity (weight) M of a given walker, is computed from the branching part 
of the Green's function Eq. 



M=G,iTZ^TZ';Ta) 



(58) 
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where Ta < t because some of the moves can be rejected. If the mean-squared distance the 
particles would diffuse in the absence of the rejection step is {r^)tot, and the actual mean- 
squared distance is (r^)^ then the actual time used in a branching step = r(r^)a/ {r'^)tot- 
In the case of multiparticle moves, the effective time step must be calculated separately from 
a computation of the accepted to attempted moves. Since M is, in general, not an integer 
one can use instead the integer 

M = int(M + , 

where ^ is a random number uniform in the range [0, 1]. In this case the average density of 
walkers is conserved and (M) = M. If M = then the walker is deleted from the ensemble; 
otherwise M — 1 copies of the configuration are made and added to the ensemble. Note that 
fixing M = 1 is equivalent to eliminate branching and, consequently, to perform a VMC 



|2 



calculation with limiting distribution f(Jl,t oo) = Hilifi'^^^l 

§ 8. After diffusing and branching all walkers in time r the mean value of the local energy 
over all walkers is computed from the obtained distribution. One has to start averaging over 
these values, after some target (equilibration) time, when the configurations are sampled 
according to the hmiting stationary distribution f(Jl,t — > oo). The target time depends 
upon the particular problem and how close $g is to the exact state. 

§ 9. Using the average value of the local energy over the whole ensemble of walkers 
{El(JV)) f(^t^oo), the trial energy Et is updated according to Ft = {Et + El) /2. Mixing this 
estimate with an old value of the trial energy, allows one to improve the convergence. 

§ 10. After computing the average energy over a sufficiently long time [rNm, where 
Nm is the number of moves per configuration), its value is stored and the first block is 
completed. Nm should be large enough for there to be little statistical correlation between 
energy subaverages obtained in different blocks. A new ensemble of iV^ walkers is generated 
by randomly copying or deleting configurations, and steps 4 to 9 are repeated completing 
another block. One has to do as many blocks as it takes in order to reach the desired 
statistical accuracy. Notice, however, that because the propagator is only accurate to C(t^), 
the distribution f{TZ,t —>■ oo) and the resulting estimates will have a time step bias. The 
way to eliminate this bias is by extrapolating all computed expectation values to r — > 0. 

Figure shows a schematic flow diagram of the algorithm presented in this Section. 



V. EXAMPLE: ELECTRON-MONOPOLE IN 



As an example application of the method we have developed in the previous Sections 
we consider the problem of a single particle of charge e, mass m* and vector position r = 
{x^, x^, x^) in IR^ confined to the surface of a sphere of radius R centered at the origin {Ai = 
S^, = 1) moving in the presence of the vector potential of a Dirac monopole at the origin. 
This problem can be solved in closed form and so constitutes an ideal model system for 
testing the accuracy of the stochastic solutions we can derive using the formalism developed 
in previous Sections. Moreover, the case of interacting electrons confined to S^ in the 
presence of a monopole field serves as the basic model which captures the essential physics 
of the quantum Hall effect ||22|| . 

The Pauli Hamiltonian for a spinless particle in S^ is 
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(59) 



where r = r/i?, and A is the monopole vector potential (V A A = B r, B being the 
strength of the radial field.) Therefore, the total number of flux quanta 2S piercing the 
surface of the sphere is given by 2S = AttR'^B /(po, where (pQ = hc/\e\ is the elementary 



flux quantum. Following Wu and Yang we can construct angular momentum operators 



L = r A {—ihV — (e/c) A) + hS f in terms of which the Hamiltonian reads 

^ IT 12 fc2o2 

If we choose a gauge where the vector potential is A = —BR cot 6 0, then the Hamiltonian, 
Eq. can be written as 



-de - dl - cote de + 2tS ^ + cot" 9 

sm 6 ^ smd 



(61) 



in terms of the usual spherical angles 6 and (y9(0<e<7r, 0<(y9< 27r, see Fig. 0.) The 
eigenstates of this Hamiltonian are monopole harmonics (normalized to 1) |2l 



>5,n,™ = Msnm (-1)^+""™ expH^^] t;^"™ ^(|w|, \v\) 

2S + n \ 



fe=0 \ I \ 




/25 + 2n + l (5 + n-m)! (5 + n + m)!y/' 
- [ n\ i2S + n)\ j ' ^^^^ 

where u = cos(e/2) exp[z(y9/2] and v = sin(e/2) exp[— z(y9/2] are spinor coordinates, n is the 
Landau level quantum number, and m = —S — n, —S — n + 1, ■ ■ ■ ,S + n is the [L^a) angular 
momentum quantum number which labels degenerate states within the n}'^ level. In the sum 
above the binomial coefficient vanish when /? > a or /? < 0. For a given S and m, the 
ground state (n = 0) and first excited state (n = 1) are 

V-gs oc u'+'^v'-'^ , (63) 
iJes OC [2{S + l)vv - (5 - m + 1)] iPg, , (64) 

respectively. The energy of a state with Landau level quantum number n is given by 

„ f nin + 1) \ TiuJr , X 

= f 2n + 1 + ^ ^ M , (65) 

where lOc is the cyclotron frequency [uj^. = \e\B /m*c). 

Let us now reformulate the electron-monopole problem in a way consistent with the no- 
tations introduced in the previous Sections. In this way one can compare the exact result 
to the numerical one obtained with the algorithm developed in this paper, thus testing the 
numerical technique. First, instead of the spherical angles 6 and (p we introduce new coordi- 
nates z and z, where z = tan(6'/2) exp[—iip], and z is its complex conjugate. Geometrically, 
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this transformation can be viewed as a stereographic projection of the sphere onto the plane, 
as illustrated in Fig. ^. The Hamiltonian can be rewritten as 

j n c 

IH7^ = — 9-'^\p. - A{z)) {p, - A{z)) g'/' + — , (66) 
in terms of the (non-hermitian) canonical momenta Pz and Pz-, and 

A{z) = -i—z ' , (67) 



2 



with metric tensor 



r{zrz)=[Z.,. ^ . (68) 

Naturally, the particles moving in the projected plane are in a space with curvature, corre- 
sponding to that of the sphere. Notice that the metric tensor is diagonal when written in 
terms of (^\ i^), such that z = + z i.e. g^''{i'', i^) = ^^^^fP^ S''" (i.e, it corresponds to 
the conformal gauge). Then, the drift is simply = D F^. 

The stochastic method developed above allows one to obtain the exact energy eigenvalues 
of the electron-monopole problem in iff we know the exact phase of the eigenf unctions. 
In other words, if the trial state is chosen such that it has the exact ground state phase, 
then independently of its modulus our stochastic approach will lead to the exact ground 
state energy. Similarly, if the trial state has a phase corresponding to an excited state 
eigenfunction then we will obtain the exact excited state energy eigenvalue. In the next two 
subsections we construct simple trial states for the ground and first excited states of the one 
particle problem. Their modulus are then used as guiding functions $g- Using these trial 
states we will apply our technique and illustrate the main ideas of our method. 

Ground State (n = 0): The 25+1 degenerate ground states of the electron-monopole 
system are labeled by their L^s = h[zdz — zdz] angular momentum quantum numbers m = 
—S, ■ ■ ■ ,S. Here we consider the m = S ground state for which the exact (unnormalized) 
wave function can be written 





\z\ 




^(1 + 


\z\ 


?) 




i',s=[ ' =|^,.|e^^-. (69) 



To illustrate our method we imagine that we do not know this exact ground state '(\}gs but 
instead we have constructed the following two trial states 



and 

= [wrm) (^)° " • 

where A and a are real valued constants. 
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The trial states i/jti and iIjt2 have been constructed so that for A = and a = they 
are both equal to the exact ground state, ipgs- For A 7^ 0, the modulus of ipxi is no longer 
equal to that of the exact ground state, but the phase of the wave function is exact, i.e. 

IV'Tll 7^ \tpgs\ , ^Tl = Vgs ■ (72) 

In contrast, for a 7^ 0, the modulus of is exact, but the phase is approximate, 



I^T2| = 1^9 



^T2 7^ Vgs 



(73) 



It follows that if ipTi is used as a trial state in a FP DMC simulation the resulting energy 
will be the exact ground state energy Eq = huJc/2, while if ipT2 is used as a trial state 
the simulation will not lead to the exact ground state energy, but instead will provide a 
variational upper bound. 

As has already been emphasized the trial state used in a FP DMC simulation should be 
constructed to be the best available approximation to the exact eigenstate, since the quality 
of the trial state can greatly influence the speed of convergence and the statistical accuracy 
of the result of the simulation. This can be clearly illustrated by considering the trial state 
iIjt2- Since the modulus of ipT2 is exact the drift velocity F which depends only on the 
modulus of the trial state will correspond to the exact drift velocity and in the absence of 
the branching term will lead to the exact density distribution. It is straightforward to show 
that 



AS 

1 + l^p 



45 e 



1 + \z\ 



(74) 



where = d^t^ In |'?/'T2p; indicating that walkers are guided away from the regions where 
the wave function is small and, in this way, the particle tends to spend most of the time 
near the top of the sphere = 0). A potential problem appears when we consider the local 
energy, 



huJr 



T2\ 



1 + 



+ b 



2\2 



aS 



S 



1 + z 



+ 



a 



4kl 



(75) 



which, of course, is not exact due to the approximate phase of the trial state. In particular. 
El diverges as \z\ —>■ 0. As we have just shown, the drift will tend to push the particle 
towards z = 0, leading to large fluctuations in the local energy. This in turn can lead to 
huge fluctuations in the population size (number of walkers), since the size of the population 
depends exponentially on the local energy. Thus, in this particular example, one has to take 
small values for a (a << 1) in order to assure fast convergence and good statistical accuracy. 

Figure ^ shows the results of FP DMC simulations, using the algorithm developed in this 
paper, for the difference between computed and exact ground state energies for trial state 
ipTi (circles) and trial state iIjt2 (squares) using different values of the time step r. The 
r — extrapolated values are also shown. For trial state ipTi we used A = 1, the number 
of walkers was chosen to be Ny^ = 300, the number of Monte Carlo steps per walker was 
2 X 10'^, and the acceptance rate was between 97% and 99.5%. For trial state 1^x2 we used 
a = 0.001, the number of walkers was A^„, = 100, the total number of Monte Carlo steps 
per walker was 10^, and the acceptance rate was the same as for ipTi- Since the parameter 
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a in '\\)t2 was chosen so that a ^ 1, the trial state and the corresponding Green's function 
were nearly exact and we were able to reach reasonable statistical accuracy with a relatively 
small number of Monte Carlo steps. As expected we find that when trial state '^ti is used 
the extrapolated energy agrees within statistical accuracy with the exact result, but when 
we use trial state i\}t2^ for which the phase is not exact, we obtain a variational upper bound 
for the exact ground state energy. 

In Fig. ^ the density profiles for the exact ground state ^ligs = ipTi{^ = 0) (Exact), the 
trial state Vti(A = 1) (VMC), the density obtained in FP DMC with trial state ^/'ti(A = 1) 
at time step r = 0.001 (FP mixed estimator), and the extrapolated density defined as ratio of 
the square of FP density to the variational density corresponding to ipTi{^ = 1); are shown. 
Note that since the density in our DMC calculation is determined as a mixed estimate (see 
Eq. |T5|), and the density operator does not commute with the Hamiltonian between the 
DMC solution and the trial state, the corresponding density profile (FP mixed estimator) 
improves on the variational result but still differs from the exact density. The extrapolated 
estimator for the density constructed by combining both, the FP mixed estimator and the 
variational density makes it possible to improve on the FP density and is seen to be very 
close to the exact result. 

First Excited State (n = 1): As a further demonstration of the validity of our method we 
turn to the first excited state of the electron-monopole system. Again, we specify the Lr^a 
angular momentum quantum number and take m = S + 1 ioi which the exact excited state 
wave function is 

/ II \S+1 
= n I l2^ 1^1 = l^-|e^^- . (76) 



We then introduce two new trial states 



5+1 

^Ti = { -7T^^ I = I^Tile^^- (77) 
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z{l + 


\z\ 
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\z\ 









and 

with the property that for A = and a = they each reduce to the exact excited state, ipes- 
As before, for A 7^ and a 7^ the modulus of V'ti is approximate and the phase is exact 
([■^Til 7^ li^esl, fTi = fes) and the modulus of ipT2 is exact and the phase is approximate 

{\lpT2\ = \lpes\, ^T2 7^ ^es)- 

Figure ^ shows the difference between FP DMC energies computed using trial states 
ipTi and il)T2 and the exact excited state energy for different values of time step r as well 
as the extrapolated r = result. The parameters (number of walkers, number of Monte 
Carlo steps, etc.) used for these simulations were the same as those used for the ground 
state simulations except that we took a = 0.0015 in ipT2- Again, our simulations gave the 
expected results — when the phase of the trial state is exact we obtain the exact energy 
(circles), Ei = (3/2 + l/S)huJc, and when the phase is approximate we obtain a variational 
upper bound on that energy (squares). 
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In Figure || the density profiles corresponding to the trial state iprii^ = 1) (VMC), 
the FP density using the same state at time step r = 0.001 (FP mixed estimator), the 
extrapolated density, computed as above by taking the ratio of square of the FP density 
and the VMC density (Extrap. estimator), and the exact density (Exact). The results are 
qualitatively similar to those for the ground state - the FP estimator improves on the VMC 
result, and the extrapolated density is nearly equal to the exact excited state density. 

The simulation results presented in this Section provide a simple test of both the FP 
DMC method and the method developed in this paper for dealing with quantum corrections 
due to curvature. The results clearly show that these methods can be used to study quantum 
systems on curved manifolds. 

VI. DISCUSSION AND CONCLUSIONS 

In this paper we have introduced a stochastic method to solve the many-body Schrodinger 
equation on curved manifolds. This method is essentially a generalized Diffusion Monte Carlo 
(DMC) technique allowing one to deal with the effects arising from the space curvature. 
We have shown that due to the curvature the diffusion matrix and drift vector, which 
appear in the Green's function used as a conditional probability in DMC simulations, acquire 
additional terms, the so-called quantum corrections. The explicit expressions for a general 
metric tensor are worked out in detail. Since the presence of the curvature leads to a number 
of other nontrivial modifications, we have presented a step by step algorithm which can be 
used to implement a code dealing with DMC simulations in curved space. 

It is worth emphasizing that our method can be applied to a wide variety of inhomo- 
geneous systems (e.g., inhomogeneous semiconductors with a position dependent effective 
mass), not just systems on curved manifolds. The reason for this, as discussed in Section III, 
is that the quantum corrections to the Green's function can be interpreted as being due to 
those terms which appear in the generalized diffusion equation describing the system once 
each of the derivatives in that equation have been commuted all the way to the left of each 
term. This definition of quantum corrections is quite general and can be applied to any 
differential equation which is second order in space and first order in time, regardless of the 
number of dimensions or any spatial inhomogeneity in the system. 

To illustrate the general methodology we have concentrated on the problem of interacting 
fermions in external electromagnetic potentials. In this case a variational upper bound to the 
exact ground state energy can be found by applying the Fixed-Phase approximation, where 
the fermionic problem is treated as a bosonic one by fixing the phase of the many-body wave 
function (which is complex-valued in general) by some trial phase. As an example, we have 
considered the problem of a single electron confined to the surface of a two-sphere, which 
has a magnetic monopole at its center. The electron thus moves in a space with curvature in 
the presence of a magnetic field which breaks time-reversal symmetry. This simple problem 
can be solved in closed form and, therefore, we have used it as a playground for testing 
our technique. In the paper we have presented two calculations, where the ground and first 
excited state energies are computed using the exact phases, but approximate modulus for 
the corresponding guiding functions. We have shown that the exact energies are reproduced 
within statistical accuracy thus proving that the approach for dealing with the quantum 
corrections is valid. 
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As emphasized in the Introduction, the method presented in this paper for performing 
DMC simulations on curved manifolds can be used to study many interesting physical sys- 
tems. An important example is the quantum Hall effect, a phenomenon which occurs when 
a two-dimensional electron system is placed in a strong magnetic field. As first pointed out 
by Haldane |^ , the electron- monopole system described in Section V provides a convenient 
geometry for performing finite size numerical studies of quantum Hall systems when many 
interacting electrons are placed on the sphere. This is in part because the spherical geom- 
etry has no boundary so that finite size effects are suppressed. In addition, the spherical 
geometry is conceptually simpler than the (flat metric) torus geometry, which also has no 
boundary, because the topological order exhibited by quantum Hall states leads to certain 
nontrivial degeneracies on the torus . 

Recently we have used the method developed in this paper to study some of the exotic 
excitations which occur in quantum Hall systems, specifically the fractionally charged quasi- 
particle excitations of the fractional quantum Hall effect PH], and the charged spin texture 



excitations (skyrmions) of the integer quantum Hall effect p6|. Previous numerical studies 



of these excitations have been based on either VMC or exact diagonalization calculations 
which, for the most part, have assumed that the wave functions describing the excitations 
are confined to the lowest (n = 0) Landau level. In fact, this is a rather poor approximation 
for real experimental systems which can exhibit significant mixing of higher Landau levels 
due to the electron-electron interaction. Because the FP DMC method allows one to go 
beyond the lowest Landau level approximation it can be used to study the effect of Landau 
level mixing on quantum Hall states and, by employing the method described in this 
paper, we were able to perform such studies using Haldane's spherical geometry. These 
simulations have provided useful quantitative results for various properties of quantum Hall 
excitations for realistic experimental parameters P5|,PB[. Along with the more rigorous test 



case of the electron-monopole system described in this paper, these calculations of Landau 
level mixing effects in quantum Hall systems using the Haldane sphere have shown that the 
method we have developed for performing FP DMC calculations on curved manifolds works 
well and can be used to study many other interesting physical systems. 
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FIG. 1. A schematic of the Fixed-Phase method for curved manifolds with general metric g^'^ . 
See the text for notation. 
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FIG. 2. Spherical and stereographic projection coordinates. R is the radius of the two-sphere 
S^. Notice that points on the sphere are projected into the complex plane {z G C) from the southern 
pole. 
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FIG. 3. The difference between computed and exact ground state energies for the trial state 
with the exact phase V'Ti (circles) and the trial state with an approximate phase ipT2 (squares) for 
various values of time step r. The r = extrapolated results are also displayed. Using a trial state 
with the exact phase in the FP DMC simulations allows one to solve the problem exactly, while 
using a trial state with an approximate phase allows one to obtain a variational upper bound for 
the exact solution. 
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FIG. 4. Ground state density for the exact ground state ipTii^ = 0) (Exact), the trial state 
ipTi{^ = 1) (VMC), the density obtained in FP diffusion Monte Carlo with trial state V'ti(A = 1) 
at time step r = 0.001 (FP mixed estimator), and the extrapolated density defined as ratio of 
the square of FP density to the variational density. The diffusion Monte Carlo density (FP mixed 
estimator) improves on the variational result but still differs from the exact one. The extrapolated 
estimator for the density constructed by combining both, the FP mixed estimator and the vari- 
ational density makes it possible to improve on FP density and is very close to the exact result. 
The density is normalized in such a way that its integral over the surface of the sphere is AttR^. 
The magnetic length is Iq = ^yhc/\e\B. 



26 



0.015 I ' 1 ' 1 ' 1 ' 1 • 1 

=!!==■= ■ ■ ■ =11= 

0.010 ■ 

• 

^ 0.005 ■ 

§ 

"J° 0.000 J t 

-0.005 I ' ' ' ' 

0.000 0.002 0.004 0.006 0.008 0.010 

T 

FIG. 5. The difference between the computed and exact excited state energies for the trial state 
with the exact phase tpxi (circles) and the trial state with an approximate phase ipT2 (squares) 
for various values of time step. The extrapolated results to time step r = are also displayed. 
For the trial state with the exact phase one finds the exact solution and for the trial state with an 
approximate phase one finds a variational upper bound for the exact solution. 
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FIG. 6. Excited state density corresponding to the trial state i^xii^ = 1) (VMC), FP density 
with the same state at time step r = 0.001 (FP mixed estimator), the extrapolated density, which is 
computed by taking the ratio of square of the FP density and the VMC density (Extrap. estimator), 
and the exact density (Exact). The FP estimator improves on the VMC result, but still differs 
from the exact density. The extrapolated density allows one to improve on the FP diffusion Monte 
Carlo result and is very close to the exact. The density is normalized in such a way that its integral 
over the surface of the sphere is 4ttR^. 
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